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FOREWORD 

The  research  described  in  this  report  was  conducted  by  the  Geophysics  Labora¬ 
tory  of  Willow  Run  Laboratories,  a  unit  of  The  University  of  Michigan's  Institute  of 
Science  and  Technology.  The  work  was  performed  as  part  of  Project  VELA 
UNIFORM,  sponsored  by  the  Advanced  Research  Projects  Agency  and  monitored  by 
the  Air  Force  Office  of  Scientific  Research  under  Contract  No.  AF  49(638)- 1759. 

The  research  period  extended  from  1  June  1966  through  30  May  1970;  the  Project 
Scientist  is  Mr.  William  J.  Best. 

The  principal  investigators  for  this  project  were  P.  Jackson,  R.  Turpenlng, 
and  D.  Willis.  This  report  was  submitted  for  publication  in  July  1970.  The  Willow 
Run  Laboratories'  report  number  is  8071-33-Fj. 
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PREFACE 


The  investigation  reported  in  this  dissertation  was  conducted  at  the 
Geophysics  Laboratory  of  the  Institute  of  Science  and  Technology,  The 
University  of  Michigan,  over  a  three-year  period.  During  this  period  the 
preliminary  results  of  this  work  have  been  described  in  presentations, 
reports,  and  a  journal  article. 

On  April  5,  1968,  a  slide  presentation  of  computer-drawn  plots  was 
made  to  the  Geophysical  Advisory  Committee  of  the  U.  S.  Air  Force  Office 
of  Scientific  Research,  Alexandria,  Virginia. 

On  April  13,  1968,  a  paper  concerning  this  work  was  presented  to 
the  Annual  Meeting  of  the  Seisnological  Society  of  America,  in  Tucson, 
Arizona  (Jackson,  1968). 

A  major  report,  including  listings  and  flow  diagrams  of  computer 
programs  was  submitted  in  September,  1968  to  the  Air  Force  Office  of 
Scientific  Research  (Willis  and  Jackson,  1968). 

On  October  5,  1969,  a  paper  was  presented  to  the  Annual  Meeting  of 
the  Eastern  Section  of  the  Seisnological  Society  of  America  in  Blacks¬ 
burg,  Virginia  (Jackson,  1969). 

A  Journal  article  describing  the  application  to  a  spherical  earth 
appeared  in  the  Bulletin  of  the  Seisnological  Society  of  America  in 
June,  1970  (Jackson,  1970). 

During  this  three-year  period  numerous  progress  reports  on  this 
investigation  have  been  submitted  to  the  Air  Force  Office  of  Scientific 
Research  and  the  American  Petroleum  Institute. 
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ABSTRACT 


Simulation  of  seismic  rays  for  a  spherical  earth  and  a  flat  earth  has  been 
achieved  in  highly  complex  models.  Travel  times  and  approximate  amplitudes  of 
seismic  waves  can  be  found  for  both  two-  and  three-dimensional  models  of  portions 
of  the  earth.  In  seismology  and  other  disciplines  ray  construction  customarily  has 
been  applied  to  simplified  geometries.  It  has  been  necessary  to  assume  that  the 
seismic  wave  velocity  distribution  of  the  earth  was  relatively  uniform  and  symmetric. 

Recently,  however,  the  earth  has  been  found  to  be  more  complex  and  non- 
uniform  than  formerly  assumed.  A  need  has  thus  arisen  in  seismology  to  test  highly 
heterogeneous  models  of  seismic  velocity  distribution.  At  the  same  time  the  devel¬ 
opment  of  the  modern  digital  computer  has  provided  a  means  of  performing  the 
necessary  ray  constructions  and  numerical  calculations. 

The  problem  of  complicated  seismic  velocity  distributions  was  therefore  investi¬ 
gated  in  terms  of  the  most  appropriate  use  of  the  digital  computer.  For  this  investi¬ 
gation  a  velocity  field  was  set  up,  and  the  propagation  computations  made  for  short 
segments  of  rays  within  this  field.  Total  travel  times  are  found  by  adding  the  travel 
times  of  connected  ray  segments.  Essentially,  the  nature  of  pro|>agation  was  dupli¬ 
cated  on  the  computer,  in  that,  at  the  location  of  each  segment  along  the  path  of  propa¬ 
gation,  the  initial  condition  and  effect  of  the  surroundings  determine  the  succeeding 
direction  of  the  following  segment. 

Both  visual  aiid  numerical  results  have  shown  that  this  simulation  method  can 
be  usefully  applied  to  investigation  of  seismic  velocity  distributions  of  portions  of 
the  earth  of  any  size  or  complexity. 
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Introduction 


Possibly  no  artificial  construction  has  been  more  useful  to  science 
than  the  ancient  concept  of  the  ray.  Not  only  in  seismology,  but  in  all 
scientific  fields  concerned  with  the  transfer  of  energy,  reasoning  with 
and  construction  of  rays  have  resulted  in  fundamental  advances. 

The  firm  basis  of  rays  in  science  is  most  dramatically  shown  in 
optics.  From  ancient  times  to  the  present,  light  rays  have  been  used 
for  intuitive  visualization  and  have  resulted  in  advances  in  optics  and 
related  fields.  Archimedes,  Aristotle,  Roger  Bacon,  Kepler,  Newton, 

Young,  Fresnel,  Rutherford,  and  Compton  all  used  the  ray  concept.  Newton's 
Optiks  (1730;  1952)  is  probably  the  most  lucid  example  of  rays  as  both 
a  rationale  and  a  means  of  intuitive  understanding  and  imaginative  dis¬ 
covery;  almost  every  figure  in  Newton's  book  is  a  ray  tracing. 

In  seismology,  with  which  this  thesis  is  concerned,  the  use  of  rays 
aided  Mohovoricic  (1910)  in  discovering  the  discontinuity  between  the 
mantle  and  the  crust;  Gutenberg  (1914)  inferred  the  existence  and  esti¬ 
mated  the  size  of  the  earth's  core  with  the  aid  of  rays;  and,  in  theor¬ 
etical  analysis,  Zoeppritz  (1919)  employed  the  ray  concept  to  derive  the 
relationships  of  reflected  and  refracted  dilatational  and  distortions 
waves  at  an  interface. 

Current  seismological  contributions  to  the  knowledge  of  the  earth 
continue  to  rely  on  the  concept  of  rays.  Current  issues  of  scientific 
journals  in  seismology  usually  contain  several  articles  which  are  il¬ 
lustrated  by  and  rely  on  seismic  rays. 


1 


In  soismics  r  ^st  ray  tracing  has  been  based  on  Herglotz-Wieche rt 
formulation  as  extended  by  Bullen  (1963).  Helbig  (1965)  used  this  geo¬ 
metrical  construction  for  a  graphical  method  for  spherical  shells. 

Sph  ically  symmetric  ray  tracing  has  also  been  described  by  Julian  and 
Anderson  (1968)  and  shown  by  Lewis  and  Meyer  (1968).  Engdahl,  et  al. 
(1968)  employed  the  formulation  for  the  1968  P-Phase  Tables  (Taggart, 
et  al .  ,  1968) . 

Iyer  and  Punton  (1963)  broke  away  from  rigid  geometrical  limita¬ 
tions  in  constructing  successive  wavefronts  by  applying  Huygens'  princi¬ 
ple  involving  complicated  logic.  Yacoub,  et  al.  (1968)  also  departed 
from  rigid  geometry  by  considering  regions  of  constant  velocity  in  which 
the  interfaces  can  be  oriented  in  any  direction.  He  also  computed  the 
amplitudes  by  Zoeppritz's  equations  for  rays  crossing  an  interface. 

The  ray  simulation  described  in  this  dissertation  is  based  upon  a 
new  concept  leading  to  the  treatment  of  heterogeneous  structures  and  the 
travel-time  solutions  for  any  discretely  specified  and/or  analytically 
represented  velocity  distributions. 

An  advance  in  the  art  of  tracing  rays  is  thus  a  significant  contri- 
oution  to  the  science  of  seismology;  and,  by  extension,  to  other  sciences. 
Currently  two  aspects  of  recent  scientific  and  technological  development 
provide,  first,  a  new  requirement  and,  second,  a  practical  means  of 
satisfying  this  requirement. 

The  requirement  to  investigate  non-simplified,  heterogeneous  velo¬ 
city  distributions,  has  arisen  from  the  discovery  that  the  earth  is  less 
synmetric  and  more  heterogeneous  than  previously  supposed.  More  accurate 
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and  extensive  seismic  measurements  are  currently  being  made.  The  need 
for  greater  understanding  has  made  simplifying  assumptions  of  the  velo¬ 
city  distribution  of  the  earth  less  useful.  Tracing  rays  through 
heterogeneous  velocity  structures  and  determining  their  travel  times 
would  be  useful  in  gaining  detailed  understanding  of  the  earth. 

The  practical  means  to  achieve  the  tracing  of  rays  and  determina¬ 
tion  of  travel  times  through  complex  velocity  distributions  is  the  em¬ 
ployment  of  the  modern  digital  computer.  Complicated  velocity  distribu¬ 
tions  imply  unpredictable  velocity  gradients.  In  large  regions  a  ray 
must  be  constructed  with  many  different  computations  because  of  changes 
within  relatively  small  regions.  With  a  digital  computer  one  can  carry 
out  successive  computations,  which  enables  one  to  treat  the  complicated 
velocity  distributions  indicated  by  modern  seismology. 

With  these  two  aspects  in  mind—the  requirement  and  the  apparent 
means  of  fulfilling  the  requirement — a  preliminary  trial  was  made  to 
develop  a  new  method  of  ray  tracing.  The  point  of  departure  was  to 
somehow  propagate  rays  through  some  type  of  mathematical  or  representa¬ 
tional  cross-section  of  a  velocity  distribution.  The  approach  was  to  use 
a  sampled  structure  into  which  a  short  segment  of  a  ray  could  be  intro¬ 
duced  at  a  given  location  and  with  which  the  directions  and  positions  of 
successively  added  segments  could  be  found  as  functions  of  the  velocity 
distribution  represented  by  the  structure. 

A  rectangular  two-dimensional  grid  of  equally-spaced  points  was 
considered.  The  grid  points  represented  positions  on  a  vertical  cross- 
section  of  the  earth,  each  point  corresponding  to  a  given  horizontal  and 
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vertical  (list  .nice  from  a  reference  point.  A  "sampled"  velocity  and  slope 
would  correspond  to  each  point  of  the  grid.  A  means  was  then  required 
ot  using  the  velocities  and  slopes,  and  the  relationships  between  neigh¬ 
boring  samples  of  their  values,  to  "propagate"  short  successive  segments 
of  the  rays  through  the  sampled  velocity  distribution.  This  means  was 
found.  An  introductory  description  follows. 

The  aforementioned  grid  was  considered  as  a  matrix  in  which  the 
rows  represented  discrete  horizontal  distances  along  the  cross-section, 
while  the  columns  represented  discrete  vertical  distances,  as  shown  in 
Figure  1.  Such  a  matrix  can  be  represented  as  a  vector  with  two  indices. 
Call  the  vector  V^. ,  in  which  i  represents  the  row,  and  J  the  column  of 
the  matrix.  For  each  set  of  indices  (i,j)  a  position  on  the  cross-section 
can  he  identified.  For  example,  let  the  interval  between  the  samples  be 
10  km.  Then  the  indices  (i  ■  1,  J  ■  1)  could  represent  the  reference 
point  (0,0)  at  the  surface  of  the  earth  and  on  the  cross-section;  the  in¬ 
dices  (i  -  10,  j  -  1)  represent  a  distance  90  km  from  the  reference  point 
on  the  surface;  the  indices  (i  ■  5,  J  ■  3)  a  horizontal  distance  40  km 
from  the  reference  point  at  a  depth  of  20  km.  The  indices  must  be  inte¬ 
gers  greater  than  zero  to  represent  the  matrix  in  the  computer.  For  this 
reason  the  smallest  pair  of  indices  (1,1)  are  taken  to  represent  an  ori¬ 
gin  (0,0).  This  bias  Is  easily  handled,  as  the  source  can  be  anywhere 
within  the  hounds  of  the  matrix,  provision  for  measuring  from  the  source 
location  must  be  made. 

At  each  location  defined  by  two  indices  there  corresponds  a  velocity 


and  a  slope 


These  values  at  discrete  points  can  be  referenced 
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from  locations  closest  to  them.  A  ray,  however.  Is  a  succession  of 
short  segments  of  possibly  differing  directions.  These  segments  may  be 
Joined  together  at  non-discrete  locations,  say  the  coordinates  (h,d), 
where  h  Is  horizontal  distance  and  d  Is  depth  in  a  cross-section  of  the 
earth.  At  each  matrix  location  (h,d)  a  computation  must  be  made  to  de¬ 
termine  the  direction  and  hence  end  position  of  a  succeeding  segment. 

The  ray  segsmnt  midpoints  could  then  be  referred  to  the  nearest  discrete  • 
(1,J)  location  to  determine  the  velocity  V  and  slope  S.  For  example, 
consider  the  scale,  or  sapling  interval  as  shown  in  Figure  1  was  10  km, 
ard  the  ray  aegswnt  midpoint  was  located  with  respect  to  the  matrix  at 
h  ■  4.75,  d  ■  3.25  corresponding  to  H  ■  37,5  km  and  D  *  22.5  km  on  the 
referenced  cross-section  of  the  earth.  The  velocity  V  and  slope  S  would 
be  found  at  the  location  indicated  by  the  indices  (1  •  5,  J  •  3),  and 
could  be  easily  referenced;  if  desired,  one  could  Interpolate  the  values 
between  those  found  at  (1  ■  5,  J  •  3)  and  (1  ■  4,  J  ■  4). 

At  the  midpoint  of  every  ray  segment  a  value  for  velocity  and  slope 

is  found.  For  geoarntrlcal  ray  tracing,  Snell's  law  requires  the  knowl¬ 
edge  of  an  incident  direction  with  respect  to  an  interface  and  the  velo¬ 
city  of  the  medium  on  either  side  of  the  Interface.  The  representational 
structure  described  above  supplies  this  information.  The  direction  of  the 
incident  ray  segment  and  the  slope  of  the  Interface  a.c  known,  therefore 

the  incident  angle  for  Snell's  law  can  be  found.  The  velocity  Vj  of  the 

medium  on  the  incident  side  of  the  Interface  is  found  as  described  In  the 
last  paragraph.  The  velocity  Vj  on  the  emerging  side  can  be  found  by  ex¬ 
tending  the  ray  segment  in  the  s«e  direction  as  the  incident  ray  for 
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the  segment,  using  the  technique  described  In  the  last  paragraph. 

For  a  fixed  rectangular  Cartesian  coordinate  system,  a  horizontal 
distance  can  he  represented  as  L  sin  0  and  a  vertical  distance  as  L  coa  0, 
where  1.  is  a  length  of  a  ray  sepnent  and  0  is  the  angle  froa  the  up¬ 
ward  vertical.  0  is  positive  in  the  counterclockwise  direction.  Now 
consider  that  the  head  of  an  Incident  ray  sagaent  of  angle  6^  located  at 
(h^.d^),  and,  by  Invoking  Snell's  law  the  merging  angle  is  found  to  be 
0  ,  as  shown  in  Figure  1.  The  location  of  the  tall  end  of  the  eaerging 
ray  segment  is  found  to  be  h^  ■  h^  ♦  0.5  L  sin  0,  d#  •  ♦  0.5  L  cos  0. 

As  the  direction  and  location  of  the  eaerging  ray  segment  as  well  as 
the  velocities  and  slopes  of  the  issMdiate  surrounding  regions  arc.  known 
whenever  they  may  be  within  the  defined  matrix,  the  eaerging  ray  sapient 
can  he  redefined  as  the  new  incident  ray  sefpent,  and  succeeding  seg¬ 
ments  generated  and  joined  end-to-end  until  a  complete  ray  is  formed. 

Travel  tines  of  the  radiation  of  seismic  waves  from  the  disturbance 
to  the  seismometer  are  of  fundamental  importance  in  seismology.  The  in¬ 
formation  to  determine  travel  times  is  available  with  this  ray  tracing 
method.  Each  segment  of  the  ray  is  located  in  a  region  in  which  the 
velocity  is  known.  The  travel  tiae  for  each  segment  can  be  found  by  di¬ 
viding  the  length  of  the  sepwnt  by  the  average  velocity  along  its 
length.  The  travel  time  for  the  rsy  is  then  the  sua  of  the  travel  times 
of  the  sepnents.  The  distance  of  travel  is  found  by  siaaslng  the  lengths 
of  the  segments.  This  distance  is  useful  for  computing  sphericsl  spread¬ 
ing  and  logarithmic  attenuation. 


6 


The  approach  described  above  was  tested  with  a  computer  program. 

The  Initial  results  were  promising,  and  the  ideas  seemed  to  be  valid. 

A  seismic  ray  could  he  simulated  through  a  matrix  representing  a  sampled 
velocity  distribution  in  a  vertical  cross-section  of  the  earth.  The  In¬ 
vestigation  to  produce  a  useful,  accurate,,  and  versatile  seismic  ray 
simulation  was  then  commenced.  The  ensuing  investigation  was  mainly  de¬ 
voted  to  attacking  the  following  problems: 

1.  Accommodating  ray  angles  through  the  entire  range  of  360*,  and 
of  using  any  defined  slope  In  association  with  any  given  ray  angle, 

2.  Reflection  from  an  interface, 

3.  Both  critical  and  non-critical  multiple  reflections, 

4.  Precisely  locating  an  Interface  between  matrix-designated  dis¬ 
crete  locations  and  extending  the  ray  to  this  Interface, 

5.  Developing  and  teating  a  method  to  obtain  high  accuracy,  of 
determining  travel  times  and  approximate  amplitude  in  a  multiply-reflec¬ 
ting  model, 

6.  Adding  the  capability  of  specifying  regions  of  velocity  charac¬ 
teristics  defined  by  an  analytical  mathematical  function, 

7.  Treating  both  curved  and  fist  earth, 

8.  Applying  the  concept  t^»  a  three-dimensional  model,  and 

9.  Making  the  method  general  so  that  one  program  can  be  used  for 
many  different  regions  of  the  earth. 

A  versatile  tool  has  been  developed  for  the  scientific  investiga¬ 
tion  of  velocity  distribution  in  the  earth.  Any  conceivable  velocity 
distribution  can  be  simulated  and  co^ared  to  actual  travel  tines  in 
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Investigations  where  geometrical  ray  approximations  are  valid.  In  addi¬ 
tion,  the  net  hod  Is  expandable,  and  can  be  uaed  aa  a  base  for  aode  con¬ 
version,  the  addition  of  Zoepprltz’s  equations  for  amplitude,  extension 
to  the  atmosphere  where  winds  provide  a  moving  medium,  and,  possibly, 
for  geometric  diffraction. 

The  simulation  of  seismic  rays  can  be  performed  with  either  two-  or 
three-dimensional  models.  Flat  earth  models  naturally  accommodate  them¬ 
selves  to  a  rectangular  Cartesian  coordinate  system.  Curved  eerth 
models,  which  at  first  glence  would  eppeer  to  be  better  treeted  in  poler 
coordinates,  are  also  treated  in  e  rectangular  Cartesian  system  because 
the  velocities  and  slopes  for  a  spherical  earth  are  easily  referenced. 
Also  the  compatibility  of  coordinate  systems  between  flat  end  spherical 
earth  enables  one  to  insert  ssmpled  values  for  anomalous  regions  within 
a  spherical  earth  model. 

As  first  conceived,  it  wee  thought  that  the  two-dimensional  techni¬ 
que  of  ray  tracing  could  be  Incorporated  into  the  three-dimensional 
technique,  and  two-dimensional  tracing  performsd  when  holding  the  values 
along  one  of  the  axes  of  the  three  distensions  constant.  However,  the 
use  of  three-dimensional  models  required  s  sufficiently  different  tech¬ 
nique  that  the  two  types  are  better  explained  separately.  The  descrip¬ 
tion  of  the  use  of  the  three-dimensional  model  is  not  self-contained. 

To  avoid  repetition  of  Identical  material,  pertinent  matter  presented  in 
the  two-dimensional  section  will  be  referred  to  when  describing  its  use 
In  three  dimensions. 
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Two-Dimensional  Simulation 


Plat  Earth 

Both  discrete  values  assigned  to  a  matrix  and  continuous  mathemati¬ 
cal  functions  may  be  employed  to  reference  velocity.  To  aid  in  visuali¬ 
sation,  the  following  exposition  is  based  on  the  discrete  case. 

Construct  a  two-dimensional  scalar  field  representing  the  seismic 
velocity  characteristics  and  slopes  along  a  plane  in  a  medium,  such  that 

V  -  f(x,s)  fl) 

S  -  g(x,s)  -  h(V)  (2) 

where  V  is  the  velocity  (actually  the  speed  with  which  seismic  P-waves 
propagate,  but  to  be  referred  to  henceforth  as  velocity,  in  keeping  with 
common  practice  in  seismology),  f,  g,  and  h  are  functions,  and  x,s  represent 
distances  parallel  to  the  axes  of  a  rectangular  Cartesian  coordinate  sys¬ 
tem.  The  functions  f  and  g  may  be  defined  differently  for  specified  regions 
of  the  field,  and  can  Include  both  continuous  mathematical  functions  and 
references  to  tables  of  arbitrary  and  discretely  changing  values.  The 
function  g  may  be  derived  directly  from  f  as  the  normal  to  the  gradient, 
or  specified  separately.  One  might  term  S  a  vector,  as  it  corresponds 
to  a  direction  and  le  based  on  the  gradient.  However,  It  is  defined  as 
an  angle  and  is  used  as  a  scalar  quantity  in  the  computation. 

Consider  a  matrix  of  discrete  values,  the  rows  of  which  represent 
equally  spaced  horlaontal  positions  within  the  scalar  field.  For  each 
element  of  the  matrix  a  velocity  and  slope  is  assigned  which  corresponds 
to  the  location  repreeented  in  the  scalar  field.  Values  for  velocities 
and  slopes  can  be  found  in  the  scalar  field  at  any  point  not  corresponding 
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to  ttu-  iii..  t.  to  locations  represented  by  the  matrix  elements.  These 
values  are  found  hv  either  using  those  of  the  matrix  elements  to  which 
the  locat ion  is  closest  or  interpolating  between  the  nearest  horizontal 
and  vertical  elements  of  the  matrix.  Thus  for  any  position  in  the  sca¬ 
lar  field,  the  value  for  the  velocity  and  slope  can  be  found  by  reference 
to  the  nearest  element  of  the  matrix. 

As  arbitrary  values  can  be  stored  in  the  matrix,  any  type  of  velo¬ 
city  distribution  can  be  referenced  to  the  scalar  field.  Through  such 
a  reference  system  a  completely  heterogeneous  velocity  distribution  can 
be  used  for  ray  tracing;  the  limits  of  heterogeneity  are  limited  only  by 
the  detail  with  which  the  matrix  elements  represent  the  scalar  field. 

That  is,  by  the  distance  represented  by  two  adjacent  elements  of  the  ma¬ 
trix.  The  matrix  with  the  velocities  and  slopes  represented  as  its  ele¬ 
ments  is  necessary  to  arbitrarily  represent  velocity  distributions,  but 
not  to  perform  ray  tracing  through  the  scalar  field. 

Select  a  position  x^,  for  the  head  end  of  an  incident  ray  seg¬ 
ment  of  length  L  and  an  angle  0^(0*<8<360*) .  The  midpoint  of  this  seg¬ 
ment  is  the  location  which  is  used  to  determine  the  incident  velocity 
in  the  manner  described  above.  Reference  6^  to  the  negative  z-axis, 
counterclockwise  positive.  The  z-axis,  which  represents  depth  in  the 
earth  cross-section,  is  positive  downward.  Extend  the  ray  segment  a 

distance  L  in  the  direction  9,. 

k 

The  initial  emerging  location  is  found  by  extending  the  ray  seg¬ 
ment  at  the  incident  angle  9k,  shown  in  Figure  1,  using  the  two  following 
equat ions : 
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x^  +  L  sin  0^ 


(3a) 


x 

P 

z 

p 


+  L  cos  0, 
k 


(3b) 


The  midpoint  of  this  extended  ray  segment  is  the  location  which  is  used 

to  determine  the  initial  emerging  velocity  V  . 

m 

Sufficient  information  (position  and  angle  of  the  ray,  and  the 
velocities  and  slopes  surrounding  the  region  of  the  ray)  is  available 
to  construct  the  ray  continuation  by  Snell's  law.  The  equation  used  for 
this  purpose  is 


9m  -  ARCSINt(Vm/Vk)  a  in  (6,-S,))  +  Sfc 


(4) 


where  0^  is  the  angle  of  emergence,  is  the  emerging  velocity,  is 
the  incident  velocity,  and  S,  (-90#<SS90°)  is  the  slope  of  the  interface 
referenced  counterclockwise  positive  from  the  x-axis.  In  case  of  con¬ 
tinuous  velocity  functions  S^  is  the  normal  to  the  direction  of  the 
velocity  gradient.  The  extension  of  this  segment  is  termed  "probing" 
in  the  sense  that  one  must  repetitively  probe  for  the  proper  emerging 
angle  0^.  For  the  first  repetition  0^  replaces  0^  in  Equation  (3) 
and  Px,  Py  are  recomputed.  A  new  0ffl  is  determined  using  the  newly-found 

V  in  Equation  (4).  Equations  (3)  and  (4)  are  repeated  in  succession, 
m 

each  time  using  the  previously  determined  value  of  0^  for  0^  in  Equation 

(3).  The  repetition  is  terminated  when  successive  values  of  0  are  suf- 

in 


ficiently  close  in  value  ( | ©m  —  where  £  is  a  predefined  small 
value) . 

This  repetition  is  necessary  because  the  location  of  the  extension 
computed  with  the  incident  angle  is  in  general  different  from  the  loca¬ 
tion  computed  with  the  emerging  angle.  If  the  ray  is  penetrating  into 
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a  region  of  increasing  velocity  the  emerging  angle  between  the  ray  segment 

and  the  normal  to  the  interface  will  be  too  large.  If  penetrating  into 

a  region  of  decreasing  velocity  this  angle  will  be  too  small.  These 

improperly  evaluated  angles  are  found  because  the  tentative  ratio  V  /V. 

m  k 

may  be  a  different  value  from  the  ratio  as  found  when  extending  the  ray 
with  the  correct  emerging  angle.  One  can  asymptotically  converge  to  the 
correct  emerging  angle  bv  repeating  these  computations.  These  repeti¬ 
tions  are  performed  to  improve  accuracy  for  computations  using  any  given 
ray  segment  length  L.  Linearity  of  the  velocity  function  is  assumed  with¬ 
in  the  distance  L.  Also,  to  improve  accuracy  one  might  reduce  the  seg¬ 
ment  length  L,  as  the  accuracy  is  an  inverse  function  of  L,  and  L  is  not 
required  to  be  equal  length  for  all  ray  segments.  However,  it  still 
must  hold  that  the  difference  between  two  successively-determined  emer¬ 
ging  angles  0  would  have  to  be  a  small  value, 
m 

When  the  sufficiently  accurate  0  is  found  by  repeatedly  applying 

m 

Equations  (3)  and  (4),  the  emerging  ray  segment  and  angle  is  determined. 
This  emerging  segment  and  angle  is  then  treated  as  the  incident  segment 
with  a  new  location  and  angle.  Segments  are  repeatedly  computed,  the 
tail  of  the  "emerging"  segment  joined  to  the  head  of  the  "incident"  seg¬ 
ment  repeatedly  until  the  ray,  which  is  the  sumnation  of  all  the  seg¬ 
ments,  penetrates  to  the  surface  or  the  boundary  of  the  defined  cross- 
section  of  the  earth. 

Spherical  Earth 

Consider  a  velocity  distribution  with  depth  for  a  symmetric  earth. 
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One  can  compute  the  velocity  and  slope  for  any  position  (x,y)  within  a 
circle  representing  the  cross-section  of  the  entire  earth.  The  velocity 
is  found  as  follows: 

Vk  -  f (R)  =  f[(xk2  +  yk2)'2]  (5) 

where  f  is  a  function  of  the  radius  R;  \?k  is  found  for  any  location  de¬ 
fined  by  (xk»yk)  from  either  a  table  and  interpolating,  or  by  computing 
some  mathematical  function  of  the  radius  R  ■  (x2  +  y2)  2. 

For  each  location  defined  by  (xk,yk)  a  slope  Sk  can  be  found, 

Sk  -  ARCTAN  (yk/xk)  ±  n/2  (6) 

where  y  and  x  are  the  coordinates  of  a  given  location,  y*0,  x*0  are  the 
coordinates  of  the  center  of  the  earth,  and  the  sign  of  tt/2  depends  upon 
the  quadrant  in  which  y  and  x  are  found. 

The  angular  distance  at  which  a  ray  emerges  is  found  by 

A  -  ARCTAN  (y^x^  (7) 

where  the  proper  quadrant  is  determined  from  the  relationships  of  the 
signs  of  x  and  y,  and  where  x  and  y  are  located  within  a  small  predefined 
distance  from  the  surface  of  the  earth. 

The  depth  of  penetration,  or,  correspondingly,  the  minimum  radius 
of  the  ray  upon  reflection  is  found  by  computing  the  minimum  radius  along 
a  given  ray,  and  retaining  that  radius  to  correspond  with  the  other  pre¬ 
viously  described  data  of  the  emerging  ray. 

Any  region  which  is  defined  by  minimum  and  maximum  values 
or  functions  of  x,  y,  or  R  can  be  postulated  so  that  a  different  function 
can  be  invoked.  For  example,  one  might  wish  to  investigate  the  travel 
times  of  rays  which  penetrate  an  anomalous  region  in  an  otherwise  symmetric 
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region  of  the  eartti.  Within  the  limits  defined,  any  sampled  value  can 
he  inserted  for  the  positions  defined  by  x  and  y,  as  in  the  flat  earth 

case. 


Three-Dimensional  Simulation 

Construct  a  three-dimensional  scalar  field  representing  the  seismic 
velocity  characteristics  of  a  medium,  such  that 

V  -  f(x,y,z)  (8) 

Sx  *  g(x,y,z)  (9) 

S2  -  h(x,y,z)  (10) 

where,  as  in  the  two-dimensional  fields  V  is  the  velocity,  represents 
the  slope  with  respect  to  the  z-axis,  Sj  the  slope  along  planes  parallel 
to  the  x,y  plane,  and  f,  g,  and  h  are  functions.  Positions  within  the 
field  are  referenced  in  a  three-dimensional  orthogonal  Cartesian  coordi¬ 
nate  system.  Directions  within  the  field  are  referenced  in  a  spherical 
coordinate  system,  in  which  the  angle  from  the  positive  z-axis  is 
4>(0<<t><Tr)  and  the  direction  from  the  a-axis  in  a'plane  parallel  to  the 
x,y  plane  is  the  angle  0(O<9<2tt)  ,  the  same  as  that  described  for  two  di¬ 
mensions.  A  segment  L  oriented  in  a  direction  defined  by  <p  and  9  pro¬ 


jects  a  distance  P  on  the  z-axis 
z 


and  on  the  x-axis  as 


and  on  the  y-axis  as 


P  ■  L  cos  <f> 
z 


P  ■  L  sin  (J>  sin  9, 
x 


P  ■  L  sin  d>  cos  9. 

y 


(ID 


(12) 


(13) 
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To  visualize  the  role  of  the  two  slopes  f-tr/2<S  j  <tt/ 2 )  and  S,(0*  S2<2ti) 

more  clearly,  consider  an  interface  or  a  plane  normal  to  the  gradient 
which  lies  obliquely  to  all  of  the  axes.  Then  the  slope  S^  is  the  di¬ 
hedral  angle  between  the  z-axis  and  the  plane  of  the  Interface 
(when  interface  is  mentioned  it  also  refers  to  the  normal  plane  to  the 
gradient).  The  slope  is  found  by  the  angle  of  the  line  which  defines 
the  intersection  between  the  x,y  plane  and  the  plane  of  the  interface, 
and  is  referenced  as  in  the  two-dimensional  case. 

The  coordinate  system  is  a  right-handed  system  in  which  the  z-axis 
is  downward.  Thus  when  a  region  of  the  earth  is  represented  positive  z 
values  correspond  to  depth  from  the  surface  of  the  earth,  or  from  an 
arbitrary  horizontal  boundary.  Although  an  orthogonal  Cartesian  coordi¬ 
nate  system  is  employed  for  position,  a  new  position  can  be  determined 
by  employing  a  length  L,  and  the  angles  0  and  6.  The  two  coordinate 
systems  are  complementary. 

Ue  see  from  the  foregoing  definitions  that,  as  in  the  case  of  the 
two-dimensional  system,  we  have  sufficient  information  to  compute  Snell's 
law  In  «"he  three-dimensional  system:  L,  0,  9,  x,  y,  z,  V,  S^,  and  S2> 
However,  in  the  three-dimensional  system  the  algebraic  and  trigonometric 
processes  become  much  more  involved. 

Consider  a  velocity  distribution  in  which  the  gradient  is  everywhere 
parallel  to  the  z-axis  as  defined  above.  Such  a  distribution  would  cor¬ 
respond  to  a  region  of  the  earth  in  which  all  interfaces  were  parallel 
to  the  surface  of  the  earth.  Therefore,  no  velocity  change  would  occur 
along  a  horizontal  direction.  With  such  a  velocity  distribution,  the 
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original  tn,-.lo  ••  of  a  rav  would  never  change,  although  its  sign  Bight. 

This  is  hot  .ms»*  V  •  V,  and,  therefore  0  ■  9.  ,  as  aeen  by  reference 
n  k  *  a  k*  7 

to  Fquat ion  (4).  The  straightness  of  the  ray  is  illustrated  in  the 
orthographic  projection  in  the  x,y  plane  of  Figure  8,  The  total  change 
in  direction  wili  be  in  the  angle  4,  which  can  be  coaputed  as 

*m  •  ARCS  IN  |  (V^)  sln(*k-Su))  ♦  Su  (14) 

where  the  subscripts  a  and  k  refer  to  eaergence  and  incidence  as  in  the 
two-dimensional  case. 

Uc  now  consider  a  scalar  field  in  which  the  interfaces  aay  be  any 
orientation  in  three-dlaensional  space;  the  condition  described  in  the 
last  paragraph  would  not  hold  for  this  case.  If  the  coordinates  were 
oriented  in  such  a  aanner  that  the  t-axis  was  perpendicular  to  the  inter¬ 
face,  then  the  condition  holds  where  6  does  not  need  to  be  recoaputed 
and  a  simple  expression  for  ♦  ,  Equation  (14),  can  be  employed.  At  any 
given  location  within  the  field  the  angles  6k  and  $k  are  known,  as  well 
as  the  a  lopes  Slk  and  S^.  If  a  coordinate  systea  can  be  found  in  which 
the  z-axis  is  perpendicular  to  the  Interface  defined  by  S^k  and  S2k*  then 
the  calculation  of  Snell's  law  for  ^  in  Equation  (14)  needs  only  to 
he  used. 

Such  a  coordinate  systea  can  be  found  by  rotating  the  axes  as  a 
function  of  SJk  and  S2k«  After  rotation  the  positions  x,y,s  are  referenced 
bv  x',y',z'  and  the  angles  $  and  6  by  and  8'  to  the  new  orientations 

of  the  coordinate  system  axes.  Snell's  law  can  then  be  invoked  with  para¬ 
meters  referenced  to  the  new  coordinate  systea.  The  new  position  is 
found  as  fellows. 
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Let  S'  be  the  angle  by  which  the  x -  and  y-axes  are  rotated  around 
the  t-axis,  and  Sj  the  angle  by  which  the  x'-axls  ia  rotated  about  the 
y'-axis  (or,  viewed  differently,  the  angle  with  which  the  x-axis  is 
tilted  to  become  normal  to  the  Interface.)  For  a  continuous  mathemati¬ 
cal  function  angle  Sj  ia 

SJ  -  (a/2)  -  arct an (AX/AY)  (15) 

where  AX  and  AY' ere  directional  derivatives  along  the  x-axis  and  v-axia 
respectively.  Sj  la  found  similarly  as 

SJ  -  (a/2)  -  arccoe(AZ/AXI  ♦  AY*  ♦  AZ1)1*  (1ft) 

where  AZ  la  the  directional  derivative  along  the  s-axls.  For  sampled 
functions  S|  and  *r* 


S'j  -  (i»/2)  -  Sjd.j.k) 

(17) 

and 

s'2  -  S2(l,J,k) 

(18) 

respectively. 

The  coordinate 

axes  can  then  be  rotated,  and  corresponding 

locations 

In  the  new  coordinate  system  found  by  rotating  first  around  the 

s-axis 

and  next  around  the 

y'-axls : 

X**  ■  X  coo  (S'2)  ♦  Y  sin  (Sp 

(19) 

Y'  •  Y  cos  (Sp  -  X  sin  (Sp 

(20) 

X'  •  X”  cos  (Sp  -  Z  sin  (Sp 

(21) 

Z*  -  Z  cos  (Sp  ♦  X”  sin  (Sp 

(22) 

The  angles  8*  and  8'  ere  found  in  a  similar  manner  by  computing  a 
position  with  trigonometric  functions,  rotating  the  coordinate  axes,  and 
taking  Inverse  trigonometric  functions  of  the  angles  of  e  directed  vector 
from  the  origin  to  the  position  In  the  new  coordinate  system. 
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Wtu-a  <.uh  segment  of  the  my  is  computed,  the  new  position  can  be 
found  bv  rotating  the  coordinates  hack  into  the  originally  specified 
coordinate  system.  From  this  position  in  the  given  coordinate  systea 
the  incident  velocity  is  known,  and  by  repeated  caaputation  in  the 
manner  described  in  the  section  on  tvo-diaenslonal  tracing  the  proper 
emerging  angles  0  and  $  can  he  asymptotically  approached  to  within  a 
predetermined  small  value. 

In  returning  the  rotated  coordinate  eyetea,  the  values  V'  and  V* 

k  a 

required  to  compute  Snell's  law  to  dataralna  tha  emerging  angle  are 

V'  -  Vk  (23) 

/•  -  Vk  ♦  (V,  -  Vk)  cos  (♦•)  (24) 

where,  again,  the  prime  refers  to  values  in  the  new  coordinate  systea, 

and  <?'  refers  to  the  repetitively  probed  value  to  correctly  determine  V  . 

m  a 

In  the  rotated  coordinate  systea  we  are  concerned  with  only  one 
direction  of  reflection,  so  this  can  ba  handled  in  the  saae  Banner  as  in 
F.quatlon  (30)  of  the  next  section.  Die  only  angle  affected  is  ♦' ,  as  8' 
will  he  unaffected.  However,  in  the  originally  given  coordinate  eyetea, 
the  direction  of  the  ray  aay  reverse  itself  In  its  projection  on  the 
x,y  plane. 


Auxiliary  Computations 

Travel  Times 

As  each  segment  of  the  ray  is  found  and  added  on  to  the  preceding 
segments,  the  travel  time  can  be  determined.  The  time  of  travel  for 
each  segnent  is  simply  the  length  of  the  individual  segaent  Lg  divided 
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by  the  average  velocity  V  at  which  It  la  propagating,  ao  thac  the  total 
travel  tlae  TT  of  a  ray  la 

n 

TT  -  l  L  /V  (25) 

a-1  *  * 


where  n  la  the  total  nuiber  of  aegaenta  froa  the  tay'a  origin  to  lta 
eaergence  on  the  aurfece. 


The  total  distance  D  a  ray  haa  travailed  la 

n 

0  -  l  L  (26) 

a-1  * 

Knowing  this  dlstanca,  one  can  ccnpute  the  amplitude  dlaenalon  due  to 
spherical  spreading  and  the  exponential  attenuation  to  energy  dissipa¬ 
tion  as 

a  As(1/D)  exp(-cD)  (27) 

where  A#  Is  the  amplitude  at  the  receiver,  A(  la  the  aaq>lltuda  at  the 
source,  and  c  la  the  attenuation  coefficient. 

In  case  of  non-crltlcally  reflected  rays  and  the  ensuing  refracted 
rays,  the  partition  of  energy  Is  rigorously  computed  by  Zoepprlts's 
equations.  The  approxlaal~ve  approach  used  here  Is  that  of  Fresnel's 
reflection  coefficient  at  nornal  Incidence.  At  a  reflective  interface 
where  non-crltlcal  reflection  occurs  the  approxlnate  reflected  aa^lltude 


A,  *  VTk  -  Vi>'<Tk *  W 

end  the  subsequently  re  free  ted  aplltude  ie 

Vl  *  \  ’  Ar 


(28) 

(29) 
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Re  f  I  «*c  t  ion 

Reflation  is  required  whenever  an  interface  of  sufficient  velocity 
difference  Is  encountered  by  the  ray,  or  whenever  the  critical  angle  la 
exceeded.  At  such  a  boundary  the  reflected  angle  0^  ie 

°R  "  *  "  ei  *  2Si  1  2nfl  (W 

where  n  is  -1,  0,  or  1,  ao  that  the  condition  0$8<2it  is  fulfilled. 

Multiple  Reflections  end  Refractions 

When  the  ray  is  both  reflected  and  refracted  at  an  Interface,  it 
splits  into  two  separate  rays.  To  aecoMOdate  this  fact,  and  to  trace 
out  both  ray  branches,  the  information  of  the  angle,  position,  total 
travel  time  and  distance  to  the  interface,  and  approximate  amplitudes 
computed  by  Fresnel*?,  reflection  coefficient  at  noneal  incidence  nuat 
he  stored  for  one  ray  branch  while  the  lnfoimatlon  for  the  other  is 
employed  for  further  rsy  tracing.  This  is  accomplished  by  continuing 
with  the  reflected  ray  until  it  emerges  or  strikes  another  Interface. 

In  emerging  at  the  surface  all  the  required  data,  such  as  distance, 
total  travel  time,  approximate  amplitude,  and  the  depth  of  penetration 
of  the  rav  is  available.  The  depth  of  penetration  is  found  by  comparing 
the  depth  at  all  segments  of  the  ray,  and  retaining  the  position  most 
distant  from  the  surface.  One  then  returns  to  the  previous  interfsce 
from  which  non-critlcal  reflection  occurred,  and,  uaing  the  information 
listed  in  the  last  paragraph  which  has  been  retained,  continues  with  the 
refracted  ray.  As  an  ensemble  of  lnfoimatlon  can  be  retained  from  pre- 
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vious  reflections,  on*  is  able  to  accomsndate  aa  many  multiple  reflec¬ 
tions  as  desired. 

The  multiple  reflections  may  be  accomodated  in  the  following  man¬ 
ner.  Consider  arrays  of  any  arbitrary  number  of  elements  for  each  para¬ 
meter  needed  to  construct  a  ray,  say 

H(i),V(i),0(l) ,TT(i) ,TL(1) ,  and  A(i),  (31) 

where  H  is  horizontal  position,  V  is  vertical  position,  0  is  angle  of 
the  first  refracted  segment,  TT  is  travel  time  to  the  Interface  under 
consideration,  TL  la  total  path  length  to  the  Interface,  and  A  is  ampli¬ 
tude;  the  index  1  in  each  value  in  (31)  represents  the  reflection  number. 
One  atarts  at  the  inception  of  a  ray  with  1*1.  At  the  first  non-critlcal 
reflection,  the  refracted  valusa  are  stored  under  the  index  1*1,  and  the 
reflected  under  the  index  i»2;  similarly,  if  another  interface  with 
non-critlcal  reflection  is  encountered,  the  ref  acted  ray  is  given  the 
index  1«2,  while  the  ref.  sy  is  given  the  :ndex  1*3,  and  so  on. 

Fbr  each  non-critlcal  reflected  ray  the  index  is  Incremented  by  1,  while 
the  corresponding  refracted  ray  is  referred  to  by  the  previous  index. 

It  le  eeen  that  any  arbitrary  niaaber  of  reflected  rays  can  be  accommo¬ 
dated  and  traced  by  incrementing  indices  in  this  manner.  After  any 
given  reflected  ray  has  reached  its  destination  the  index  i  is  reduced 
by  1,  and  the  previously  refracted  ray  is  then  continued  by  uaing  the 
initial  conditions  of  the  parameters  listed  in  (31)  with  the  next  lower 
index  maber.  It  can  be  seen  that  all  refracted  rays  will  be  traced  out 
until  the  index  reverts  *o  1*1,  which  is  the  index  nimber  of  the  inci¬ 
dent  ray.  Also,  because  previously  refracted  rays  may  subsequently 
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encounter  .1  now  set  of  multiple  reflections,  all  multiple  reflections 
of  the  originally  reflected  or  refracted  rays  will  be  accommodated. 

Interface  l-ocatlons 

Because  interface  locations  are  so  critical  In  affecting  the  loca- 
tlon  of  the  emergence  of  a  ray,  a  means  was  developed  of  precisely  lo¬ 
cating  the  interface  and  extending  or  retracting  the  ray  to  terminate 
precisely  on  the  interface.  This  Is  accomplished  In  the  lollowlng  man¬ 
ner.  The  location  of  a  point  Is  found  In  the  square  defined  by  four 
matrix  elements  and  located  on  the  interface.  The  horizontal  and  verti¬ 
cal  distance  from  the  ray  segment  head  end  to  this  point  Is  determined. 
Knowing  the  slope  of  the  Interface,  the  angle  of  the  ray,  and  the 
horizontal  and  vertical  distances  to  the  defined  point,  a  series  of 
computations  In  which  the  law  of  sines  Is  Invoked  Is  then  made.  Many 
different  angular  conditions  are  Involved  as  a  result  of  combinations  of 
ray  angle,  slope,  and  positive  or  negative  horizontal  distances.  For 
details,  one  can  consult  the  subroutine  GINMAD  In  the  Appendix,  where 
the  computations  are  given  In  Fortran. 

Computer  Programs 

Three  main  programs  have  been  developed.  The  first  Is  for  two- 
dimensional  flat  earth,  the  second  for  two-dimensional  spherical  earth, 
and  the  third  for  three  dimensions.  All  have  been  written  In  the  most 
general  computer  language,  Fortran.  The  version  Is  Fortran  IV-G,  each 
program  of  which  has  been  compiled  and  run  on  the  IBM  360/67  computer 
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of  The  University  of  Michigan.  The  programs  include  plotting  instructions 
which  have  been  used  for  the  illustrations. 

Subroutines  have  been  developed  to  precisely  reflect  at  an  interface, 
to  extend  the  ray  precisely  to  the  surface,  to  compute  Snell's  law  for 
any  incident  angle  across  an  interface  of  any  slope,  to  increment  initial 
horizontal  and  vertical  positions  and  angles,  and  to  compute  reference 
travel  times  for  each  ray  making  up  an  onset  ting  plane  wave. 

The  indexing  for  the  spherically  symmetric  case  has  been  made  com¬ 
patible  with  that  for  velocity  distributions  with  a  rectangular  reference 
system,  in  such  a  manner  that  the  data  for  horizontal  layering  structures 
can  be  read  in  by  a  Fortran  Namelist.  In  this  way  anomalous  regions 
which  are  very  heterogeneous  can  be  Included. 

The  programs  have  been  made  general,  so  that  a  different  number  of 
samples  and  size  of  sampling  interval,  number  of  multiple  reflections. 
Incremented  or  changed  output,  number  of  Initializations,  etc.,  can  be 
entered  without  recompiling  the  programs.  Thirty  to  forty  values  are 
read  from  Fortran  Namelists;  the  values  can  be  changed  individually. 

Fach  ray  increment  involves  calling  trigonometric  functions  approxi¬ 
mately  7  to  25  times  (usually  about  15  times).  For  an  indication  of  the 
time  required,  see  Figure  7.  The  computation  for  Figure  7,  including 
loading,  20  travel-time  printouts,  and  plotting  instructions,  required 
10  seconds  of  IBM  360/67  computer  time.  It  is  felt  that  this  is  reason¬ 
able  and  economical  time  for  such  ray  tracing. 

The  three  main  programs  and  their  subroutines  are  listed  in  the 
Appendix. 
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Results 


Accuracy 

The  method  of  ray  simulation  was  developed  primarily  with  synthetic 
data.  First  very  simple  data  were  used  for  testing  and  developmental 
purposes.  As  the  work  progressed,  it  was  felt  that  a  test  was  necessary 
to  determine  its  accuracy. 

For  this  reason  the  method  for  tracing  rays  in  spherically  symmetric 
earth  was  attempted.  For  at  least  60  years  the  earth  has  been  investigated 
as  a  spherically  symmetric  structure  in  which  the  velocity  gradient  every¬ 
where  points  directly  at  the  center  of  the  earth.  Much  of  seismology 
has  been  concerned  with  this  model.  Because  of  the  continued  Investiga¬ 
tion  of  this  model  over  most  of  the  world  for  such  a  long  period  of  time, 
many  travel  times  from  Individual  disturbances  at  many  distances  have 
been  recorded. 

From  these  travel  times  a  velocity-depth  distribution  of  the  earth 
has  been  inferred.  This  distribution  is  one  of  the  fundamental  problems 
of  seismology,  and  most  seismologists  have  been  concerned  with  it. 

Recently  Herrin,  et  al.  (1968)  published  the  1968  Seismological 
Tables  for  P  Phases  in  a  special  issue  of  the  Bulletin  of  the  Seismological 
Society  of  America  (BSSA).  These  tables  represented  the  latest  refine¬ 
ments  of  the  contributions  of  the  science  of  seismology.  In  this 
special  issue  a  velocity  vs.  depth  distribution  of  the  mantle  was  given. 

As  both  the  velocity  distribution  and  travel  times  were  given  in  this 
special  issue  of  BSSA,  a  model  was  available  to  test  this  method  of 
ray  tracing. 
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As  one  might  expect,  the  ray  tracing  method  was  initially  shown  to 
be  inaccurate  in  the  first  tests.  Two  improvements  were  made  in  the 
technique  as  it  then  existed.  Tho  first  was  to  take  the  incident  velo¬ 
city  at  the  midponnt  of  the  incident  ray  segment  and  the  emergent  velo¬ 
city  at  the  midpoint  of  the  emergent  ray  segment.  The  second  change 
was  to  repeat  the  computation  of  Snell's  law,  as  explained  in  the  two- 
dimensional  modelling  section,  until  the  proper  emerging  velocity  was 
found.  It  was  found  with  a  testing  program  that  the  proper  value  of 
emerging  velocity  was  approached  in  a  damped,  oscillating  fashion.  It 
was  then  a  simple  expedient  to  compare  successive  values  until  the  dif¬ 
ference  between  them  was  less  than  an  arbitrary  value. 

This  latter  change,  for  which  the  need  was  not  formerly  apprehended 
was  the  key  to  achieving  as  much  accuracy  as  desired,  without  reducing 
the  length  of  the  ray  segments  to  require  an  uneconomically  large  amount 
of  computation.  As  Snell's  law  is  the  basis  of  both  this  method  and 
that  used  for  the  1968  P-Phase  Tables  (Engdahl,  et  al.  1968)  the  two 
methods  should  result  in  Identical  relationships  between  velocity  distri 
bution  and  travel  times.  Differences  should  be  due  only  to  the  size  of 
the  sampling  intervals,  length  of  the  incremented  segments,  and  method 
of  interpolation. 

The  equivalence  between  the  two  methods  has  been  found.  As  the 
sampling  interval  and  segment  length  were  reduced,  the  travel  times 
shown  in  the  1968  P-Phase  Tables  have  been  approached  more  and  more 
closely  when  using  the  same  velocity  distribution. 
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The  accuracy  is  illustrated  in  Figure  2,  where  the  velocity  distri¬ 
bution  used  for  the  1968  P-Phase  Tables  was  approximated  with  samples 
taken  at  15  km  intervals  with  9  km  segments.  Figure  2  is  a  reproduction 
of  a  computer  output  using  this  ray  tracing  method  and  also  shows  cor¬ 
responding  distances  and  travel  times  from  the  1968  P-Phase  Tables.  The 
average  difference  between  the  two  methods  in  rravel  times  for  the  five 
P-rays  shown  in  Figure  1  is  .252  seconds,  or  an  average  difference  in 
travel  times  between  methods  of  1/100  of  1%.  As  the  standard  deviation 
of  the  P-Phase  Tables  is  1  second,  .252  seconds  is  considered  satisfac¬ 
tory.  Every  indication  is  that  the  travel  time  vs.  velocity  distribution 
values  could  be  made  closer  by  further  reduction  of  the  sampling  interval 
and  segment  length. 

Plane  Waves 

Figures  3  and  4  were  drawn  to  investigate  the  simulation  of  a  plane 
wave.  The  velocity  structure  is  based  upon  a  hypothetical  velocity  dis¬ 
tribution  within  a  cross  section  from  the  Pacific  Ocean  across  the  Coast 
Ranges,  Central  Valley  of  California,  the  Sierras,  and  part  of  Nevada. 
Figure  3  represents  a  PcP  wave  from  A  »  30°  and  Figure  4  a  P  wave  from 
the  same  distance.  The  simple  subroutine  Planit  was  used  as  a  reference 
for  travel  time  across  the  wavefront.  These  figures  were  obtained  in 
May  1969  as  an  aid  in  Charles  G.  Bufe's  investigation  of  PcP  waves  (1969). 
Bufe  was  able  to  estimate  PcP-P  travel-time  differences,  and  anticipate 
arrival  anomalies. 
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Spherical  Earth 


Figures  5  through  8  illustrate  two-dimensional  seismic  ray  tracings 
of  an  entire  cross-section  through  the  earth.  In  Figure  5  only  a  single 
reflection  was  allowed,  enabling  one  to  obtain  PcP,  PKKP,  PKIKKIKP,  and 
P  rays.  Had  no  reflection  been  allowed  below  the  critical  angle  of  in¬ 
cidence,  only  the  PKP,  PKIKP,  and  P  would  have  been  drawn  in  Figure  5. 
With  four  reflections  allowed  one  can  generate  PKP  type  waves  up  to 
PKKKKP,  as  shown  in  Figure  6. 

In  Figure  7  a  source  at  700  km  depth  was  simulated.  Rays  were 
traced  from  this  source  at  18°  intervals,  and  each  ray  traced  inde¬ 
pendently  of  the  others.  Note  that  the  second  K  leg  of  the  initial 
18°  ray  is  overtraced  by  the  second  K  leg  of  its  symmetric  counter¬ 
part  the  360°  -  18°,  or  342°,  ray.  This  overtrace  indicates  high 
positional  accuracy,  in  that  A  and  travel  times  differ  by  less  than 
.01%  between  symmetric  rays. 

Three  Dimensional 

Illustrations  of  computer  plots  in  three  dimensions  are  shown  in 
Figures  8  through  11.  To  clearly  show  the  paths  of  the  rays  in  three 
dimensions  a  computer  program  was  devised  to  show  a  perspective  and 
three  orthogonal  projections.  Figures  8  through  11  show  a  single-point 
perspective  in  the  upper  left  portion  of  the  figures,  and  three  or¬ 
thogonal  projections  along  the  principle  planes  in  the  remaining  portions 
of  the  figures.  Consider  that  the  (x,y)  plane  represents  the  surface  of 
the  earth,  and  the  positive  z-axis  corresponds  to  depth  in  the  earth. 
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In  Figure  8  and  Figure  9  a  velocity  distribution  increasing  with 
depth  and  having  zero  slope  everywhere  with  respect  to  the  surface  is 
shown.  In  Figure  8  cones  of  rays  from  a  point  source  at  the  surface 
are  shown;  ravs  of  a  sufficiently  large  initial  angle  undergo  critical 
reflection  and  return  to  the  surface  (x,y  plane).  In  Figure  9  a  line 
source  at  depth  is  shown.  Points  are  chosen  along  the  line,  and  fans 
of  rays  are  propagated  toward  the  surface.  It  should  be  pointed  out 
here  that  any  type  of  initializing  rays  may  be  chosen.  For  example,  in 
Figure  9  the  points  of  origin  along  the  line  could  have  been  made  much 
closer  together,  and  the  angular  extent  of  the  fan  decreased,  with  any 
choice  of  angles  between  individual  rays  chosen.  As  described  in  the 
section  on  the  computer  programs,  the  values  for  the  increments  of  these 
points  and  angles  can  be  entered  in  the  Nanelist. 

Figures  10  and  11  illustrate  an  oblique  gradient  with  respect  to 
the  z-axis.  The  mathematical  distribution  of  the  velocities  in  Figure 

10  is 

V  -  6.0  +  X  +  2Y  +  3Z  02) 

These  illustrative  figures  make  use  of  velocity  values  which  exceed  nor¬ 
mal  earth  velocities.  These  velocities  were  used  for  the  purposes  of 
demonstrating  the  capability  of  tracing  such  a  distribution.  In  Figure 

11  a  sampled  layer  extending  from  a  level  corresponding  to  40  km  depth  to 
190  km  depth  is  placed  within  the  field  described  for  Figure  9.  The 
gradient  is  parallel  to  the  z-axis  within  this  slab,  and  from  140  km  to 
190  km  depth  the  velocities  are  constant.  As  shown  in  Figure  10,  one  can 
define  regions  in  three-dimensional  space  which  can  be  either  analytical 
mathematical  functions  or  sampled  data.  The  ntmber,  sizes,  or  shapes  of 
the  regions  are  not  restricted. 
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Potential  Application  to  the  Seafloor  Spreading  Problem 


In  the  theory  of  seafloor  spreading,  the  earthquake  zones  are 
associated  with  upward  veiling  zones,  such  as  the  Mid-Atlantic  Ridge, 
where  surface  material  is  being  produced,  and  the  boundary  between 
large  lithospheric  plates  as  in  the  Circum-Paclf lc  zone. 

The  problem  of  these  boundaries,  one  of  which  is  impinging  on  the 
other,  is  treated  by  Iaacks,  Oliver,  and  Sykes  (1968).  They  postulate 
that  one  plate  underthruata  the  other,  particularly  in  island  arc  re¬ 
gions.  At  this  boundary  the  lithosphere  (the  upper  100  km  of  the  earth's 
surface)  of  the  oceanic  plate  is  bent  downward  into  the  underlying 
asthenosphere  of  the  continental  plate.  The  precise  geometry  in  which 
the  two  plates  Join  is  not  known,  and  differs  between  regions.  Isacks, 
Oliver  and  Sykes  (1968)  postulated  several  configurations  as  a  function 
of  local  conditions  and  rate  of  underthrusting.  Mltronovas,  Isacks,  and 
Seeber  (1969)  have  investigated  this  problem  in  the  Tonga  Islands  arc. 
Murdock  (1969)  has  discussed  velocity  distribution  under  the  Aleutian 
Region.  Carder,  et  al.  (1967)  have  docuamnted  travel-time  anomalies 
from  the  LONGS HOT  explosion. 

Figures  12  through  15  were  drawn  with  an  arbitrarily-chosen  geometry 
based  on  the  models  postulated.  In  the  geometry  chosen  the  6.75  km/sec 
layer  is  taken  as  the  topmost  layer  which  bends  downward  and  the  surficial 
6.00  km/sec  layer  is  unaffected  by  bending.  These  figures  show  different 
earthquake  locations  and  the  directions  of  seismic  waves  propagating  from 
them.  Two  of  the  earthquakes  are  directly  above  and  two  directly  beneath 
the  upper  boundary  of  the  under  thrusting  lithosphere.  The  two  earthquakes 
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.•bow  1 1:«  '  M.:ul.ir\  are  located  in  rontons  where  volcanic  activity  and 
-uiKM  m  iiMi  u  t  ivitv  is  indicated  (Mitronovas,  lsacks,  and  Seeber, 

II.  V  t  ■„.>  i  »t  t  lujii.ikes  beneath  tin*  boundary  are  located  where  most 
seisnii  st  ii'.is  and  the  highest  probability  of  earthquakea  are  thought 
to  occur.  The  rays  for  these  plots  were  generated  at  4*  radial  incre¬ 
ments  from  a  point  source. 

Although  the  frequency  of  earthquakes  with  sources  above  is  less 
than  the  frequency  of  those  within  the  underthrusting  lithosphere,  both 
source  regions  have  been  shown  for  contrast.  This  contrast  arises  fro* 
trapped  waves  within  the  low  velocity  layer  when  the  source  is  within 
the  layer.  Rays  from  within  the  layer  critically  reflect  at  tha  interface 
at  angles  greater  than  55*  fron  the  normal,  to  the  Interface.  Raya  from 
sources  outside  the  low  velocity  layer  can  be  critically  reflected  only 
when  the  inclination  of  the  Interface  through  which  tha  raya  have  entered 
is  oriented  differently  fron  the  interface  through  which  they  would  have 
otherwise  emerged.  Parallel  Interfaces  hounding  a  low  velocity  layer 
cannot  trap  rays  whlcu  enter  fron  outside  the  layer.  In  the  model  used 
for  Figures  12  through  15  critical  reflection  fron  externally  generated 
rays  can  only  occur  where  the  interface  of  the  low-velocity  layer  is 
curved,  as  shown  by  the  one  critically  reflected  ray  in  Figure  12. 

The  effects  of  the  low-velocity  layer  are  clearly  naan  In  Figures 
12  through  15.  Shadow  rones  are  found  at  the  surface  when  the  source  is 
above  the  layer.  These  occur  because  the  rays  penetrating  the  layer  ere 
bent  toward  the  normal  of  the  Interface.  These  rays  are  displaced  from 
the  paths  they  wouldhave  taken  in  the  absence  of  the  low-velocity  layer. 


30 


and  they  emerge  at  a  greater  distance  from  the  source  than  would  other¬ 
wise  be  the  case.  From  sources  within  the  low-velocity  zone,  the  influ¬ 
ence  of  critical  reflections  is  shown  in  Figures  14  and  15.  Three 
separate  groups  of  rays  can  be  seen  propagating  up  the  low-velocity  layer: 
those  initially  reflecting  from  the  upper  interface,  those  which  do  not 
touch  either  the  upper  or  lower  interface,  and  those  initially  reflecting 
from  the  lower  Interface. 

Travel  times  from  each  of  the  four  sources  are  shown  in  Figures  16 
through  19.  In  Figures  16  and  17,  the  focal  depths  are  SO  km  and  100  km 
respectively,  with  sources  above  the  layer.  The  shadow  zones  caused  by 
the  ray  dlsplaceswnts  are  evident.  Travel- time  anomalies  above  2  seconds 
occur  near  the  lladt  of  the  shadow  cones  nearest  the  source. 

Travel  times  for  sources  within  the  downward  bending  layer  are  shown 
in  Figures  18  and  19,  where  the  focal  depths  are  60  and  110  km  respective¬ 
ly.  Critical  reflections  occur  with  a  large  angular  extent  from  the 
source.  The  three  separate  arrival  sequences  are  produced  by  the  three 
groups  of  rays  described  above.  The  three  separate  legs  of  the  travel- 
tlme  curves  occur  at  distances  approximately  the  same  as  those  where  the 
shadow  cones  are  produced  for  the  sources  above  the  layer  as  shown  in 
Figures  17  and  18. 

These  results  Indicate  that  large  effects  upon  travel  times  are 
found  within  50  km  of  the  point  where  the  under thrusting  lithosphere 
starts  to  bend  downward  from  the  unperturbed  oceanic  plate.  A  pertinent 
investigation  of  seismic  propagation  in  the  Tonga  Islands  arc 
(Mlctonovas,  et  al.  1969)  was  conducted  with  five  seismometer  locations 
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P  iced  l! 


•  I  es  ill 


!.’*»  » n  in  '00  kn  from  this  point  toward  the  under- 


t li runt  m,;  li'  m  I'luMr,  Hi«  nearest  s«'i sraone ter  location  toward  the  oceanic 
plate  w.i-.  ippr**xiia.»teiv  JJO  km  dist.int  on  Raratonga  in  the  G>ok  Islands. 

I'n  to  rt  unai  e 1  v  ,  the  loc.it  ton  of  the  sci  smiroeter  s  were  dictated  by  the 
positions  it  the  available  Isl.uids,  and  were  not  within  the  region  in 
whith  r  tv  ti  icing  stmws  shadow  /ones  and  multiple  arrivals. 

To  investigate  the  travel  times  beyond  the  region  shown  in  Figures 
1J  through  11,  the  spherical  earth  program  was  attached  as  a  subroutine 
to  the  flat  earth  program.  Many  rays  exceeded  the  boundaries  shown  in 
Figures  !’  through  11  without  reaching  the  surface.  These  rays  were 
used  as  sources  for  the  spherical  earth  program.  This  was  accompl lshsd 
hy  saving  the  location,  angle,  travel  time,  and  distance  travelled  as  an 
input  to  the  spherical  earth  subroutine.  Travel  times  out  to  approxi¬ 
mately  60*  were  found.  These  closely  paralleled  the  1968  P-Phase  Tables 
for  corresponding  depths  of  focus.  The  anomalies  found  were  a  max lata 
of  .  *i  seconds,  uid  were  at  the  maximum  in  the  region  near  25*.  No 
cone  Ills  ions  were  drawn  fieri  these  results. 

Figures  ;  J  through  1‘‘  illustrate  the  diagnostic  capability  of  the 
method  presented  if  this  study.  The  travel  times  are  particularly  re¬ 
vealing  it  distances  within  50  km  of  the  point  at  which  the  lithosphere 
starts  to  herd  downward.  They  depend  upon  the  location  of  the  source, 
vary  with  direction,  and  produce  distinctive  features  in  the  ray  dia¬ 
grams  and  travel-time  curves.  Civcn  a  set  of  seismograms  at  suitable 
distances  from  the  earthquake  in  a  continental  plate  boundary  area,  an 
invest f gat  ion  of  the  velocity  distribution  and  source  depth  can  be  aided 
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with  this  method.  A  more  detailed  investigation  would  also  Include  azi¬ 
muthal  effects,  which  could  he  found  by  projecting  the  given  model  on 
croas-sectlons  at  various  angles  with  the  cross-section  shown. 

It  is  not  the  purpose  of  this  study  to  investigate  a  portion  of  the 
earth  in  detail,  but  rather  to  present  a  method  to  aid  in  such  investi¬ 
gations.  For  example  the  travel  times  and  anomalies  found  with  curved 
earth  out  to  60*  were  not  pursued  as  to  their  meaning,  but  rather  a  method 
to  obtain  them  using  a  spherical  earth  subroutine  was  demonstrated.  Also, 
the  modal  chosen  for  the  underthrusting  lithosphere  did  not  Include  the 
6.00  km/sec  layer.  Many  variations  of  the  model  would  be  made  during  a 
detailed  investigation— variations  such  as  thicknesses  of  the  layers, 
effects  of  depth,  and,  therefore,  heat,  on  the  velocity  of  the  layers, 
angle  of  downward  bending,  extent  of  the  downward  thrusting  layer, 
effects  upon  the  overlying  layers,  departures  from  horizontal  layering, 
and  many  other  variations.  It  is  apparent  that  such  changes  can  be 
accommodated  by  this  ray  simulation  method.  The  results  of  these  changes 
are  given  both  visually  and  numerically. 

Conclusions 

This  study  has  resulted  in  an  ecortoadcal  and  versatile  method  to 
aid  in  seismic  analysis.  Geometric  ray  tracing  can  be  performed  with 
models  representing  any  conceivable  velocity  distribution  within  the 
earth.  Refractive  as  well  as  reflective  computations  can  be  performed, 
including  multiple  reflections.  Travel  times  and  approximate  amplitudes 
can  be  determined  for  these  velocity  distributions. 
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,>!  (  .  primary  goals  of  seismology  Is  to  model  the  velocity 
distribution*,  witiiin  the  earth.  Hypotheses  of  velocity  distributions 
are  »  01  tmuouslv  being  produced  by  seismologists.  An  important  test 
tor  the  validity  of  these  models  has  been  developed. 

Uunpiiter  drawn  plots  and  numerical  output  of  travel  times  have 
been  used  to  illustrate  the  capabilities  of  this  simulation  method. 

The  travel  times  can  be  made  arbitrarily  close  to  those  published  by 
t he  Se ismo logical  Society  of  America.  This  method  of  geometrical  ray 
tracing  performs  as  accurately  and  with  more  versatility  than  previous¬ 
ly  published  ray  tracing  techniques. 

It  is  hoped  and  expected  that  this  technique  will  be  used  by  the 
svismo logical  profession  and  that  it  will  be  extended  and  Improved,  aa 
anv  useful,  new  contribution  should. 
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FIGURE  1.  REFERENCE  FRAME  FOR  INCREMENTING  RAYS.  Seismic  velocity  and 
slope  are  associated  with  each  matrix  element.  Referencing  for  the  slope  is  shown  at 
element  (2,  1).  An  Interface  between  two  velocity  layers  can  be  positioned  between 
matrix  elements  as  shown  by  the  dotted  line  through  B.  The  scale,  or  sampling  inter¬ 
val,  is  the  distance  between  matrix  elements  as  represented  in  kilometers.  An  incident 
ray  segment  of  length  L  and  initial  angle  \  extends  from  A  to  B.  Coordinates  of  B  are 
found  by  extension  from  A:  Bx  =  Ax  ♦  L  sin  9y,  Bz  *  Az  ♦  L  cos  The  segment  BPj 
is  the  first  "probing"  segment  to  obtain  initial  values  of  seismic  velocity  at  midpoint  1. 
The  angle  is  computed  with  Snell's  law,  using  emerging  velocity  at  midpoint  1.  in¬ 
cident  velocity  at  midpoint  k,  and  incident  angle  "Probing"  segments  are  recom¬ 
puted  using  successive  values  of  '<m  until  the  absolute  value  of  (  >m  -  is  less  than 

a  predetermined  small  value.  The  velocity  and  slope  associated  with  the  segment  BPm 
and  the  last  computed  value  of  are  then  used  as  the  incident  ray  parameters,  and  a 
subsequent  emerging  angle  computed  for  the  ray  segment  proceeding  from  point  Pm. 
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COMPUTER  OUTPUT  |C  OR  RES  PON  DING  TRAVEL  TIMES 
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FIGURE  S.  PLOTS  OF  P,  PcP,  PKP,  PKIKP,  PKKP,  AND  PKIKK3KP  RAYS  GENERATED 
WITH  A  SPHERICAL  EARTH  PROGRAM.  Jeffreys -Bullen  core  and  1968  P-phase  mantle 
velocity  distributions.  Radial  velocity  distribution  sampled  at  70  km  Intervals.  Rays  seg¬ 
ment  lengths  70  km. 
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FIGURE  6.  MULTIPLE  REFLECTIONS.  PcP  and  PKP  rays  up  to  PKKKKKP.  Four 
multiple  reflections  preset  into  computer  program.  Any  number  of  multiple  reflec¬ 
tions  can  be  preset  into  program. 
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FliiFRE  8  THREE  DIMENSIONAL  RAY  TRACING  IN  A  REGION  DEFINED  PY  30' 
30  30  VELOCITY  SAMPLES  Constant  velocity  in  planes  parallel  to  the  x.  y 

plane,  me reas mu  velocity  in  z-direction  for  IS  sample  points  from  the  x.  y  plane. 
Four  cones  of  rays  originating  from  point  in  the  x,  y  plane.  Rays  emerging 
at  the  \,  v  plane  shown  by  slash  marks  with  orientation  ol  angle  of  emergence. 
Fpper  left  one-point  perspective;  remaining  views:  labelled  orthographic  pro¬ 
jections. 
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FIGURE  9.  LINE  SOURCE  AT  DEPTH  FOR  THREE-DIMENSIONAL  REPRESENTA¬ 
TION  AS  SHOWN  IN  FIGURE  8.  Plotting  defects  in  upper  right  (x,  z)  projection;  all 

lines  should  extend  to  x-axts. 


FIG1  RF  10  rilHKK  -dimknsional  ray  tracing  with  velocity  distribu 

MON  I)F  FINED  PY  CONTINUOUS  MATHEMATICAL  FUNCTION.  Velocity  4.0 
0  I  »\  •  o.3v  •  0.6/..  where  volume  is  defined  as  in  Figure  8. 
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Unaffected  by  Lithosphere  Bending 
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FIGURE  16.  EARTHQUAKE  ABOVE  UNDERTHRUSTING  LITHOSPHERE  AT  CON¬ 
TINENTAL  PLATE  BOUNDARY.  Depth  of  focus:  50  km.  Above:  Travel  times  to¬ 
ward  and  away  from  underthrusting  lithosphere.  Below:  Differences  between  travel 
times  toward  and  away  from  underthrusting  lithosphere.  See  Figure  12  for  ray  paths 

and  velocity  distribution. 
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APPENDIX 

Listings  of  Computer  Programs 
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»  AwPAV  lip  SAMP|  FI)  WFIliriTV  WAIOFS 


Ml||.  MtM[W|iM  AN|.|  |.  AT  WW|(  M  nPFf|t,F  |  N  I  F  W  F  A  P  F  |llPAI|OI>i  |S  IlSFM 

FI  APBAV  f  (lFlTAtMfMI-  SAMI>|  Ml  SI  OFF  VAIIIFS 

i*cr.tM  imitiai  amc.i  f  fit-  pay 

hpim  |  r*r  pfmfmtation  anc.i  f  iif  |  at  t  tif  mt  pavs 

|IFC  FXPIINFM1  I  A|  OHPFMFMt  MV  AT  TF  Ml  I A  T  1 1 IM 

PF  I  M  Ml  IP  I  /HUT  Al  tMrpFMFMTAT  1 1  »M  t)F  |  MP,  IIIFMT  PAYS  _  _  _ _ 

l»F|  V  WFPTtlAI  I  Mf  PFMFMTATI  f'M  IIF  |N(|llFMT  PAW 
IMFBAC  PISTAMCF  IN  KM  H  F  T  U  F  F  N  SAMPIF  PII|MTS 
P  PAY  SFPMFMT  I  FNC.TH 

it  vn  F  maximum  mmmmfp  of  bay  jmcbfmfnta  t I iims 

ifxact  *i«-*  i  tt  m  tii  f  Mr.i  iiiVf  pbpcT  sf  imtfbfapf'  i.  nriT'T  i  in  s 

t  MIU/ 1  F  V  I  TPM  TO  Mill  TIPI  Y  r,  MY  V  F I  Of.  I  T  Y  TO  f;  I V/ F  APPF  Ah  ANT.  F  IIF  PAY*; 


-  ...ATI  ™M  PPIiPIIP  T  1 1  IN  At  TII  UFI.riflTY  WM|-N  MAKIIvl,  A  (.1  INPUT  F  tt  MI|V 

in  amp  sw  itch  tii  |mpi  imp  mamfi.  1ST  impiit  i  if  vFinr.lTv  ami  slopf  Iiata 

ipianf  ‘-hIttm  in  ifir.i.imf  r  At  pm.  At  iom'S  fop  impToW  pTami?  wav'f 

IPIOT  ^  w  I  T  r  W  TO  t  NCI.  MOP  CAIPOMP  P|.  FIT  T  I  NP 

I  PHO<;h  I  I M I  T  OF  F'llMMPP  IIP  HOPI/ONTAt  S/.MPI  F  Pfl  I  MTS 
.IPOPFM  I  I M  t  T  OP  I/FPT!rA|  NIIMHFP  OF  VFMT  I  f  Al  SAMPIF  POINT*; 


i  m  ”  '  ' '  I  ***1 .  i,  mi  »r  ii|  nMNi'^i  ii,  ii  riiMM 

ma  lw  switch  to  imclmof  C.ont  iNoniis  matmfmatifal  function 
NMIII.T  MIIMPFP  nr  MIIITTPI'F  BPFCFr.TTfW*  mtoftpit 

NOB  NIIMMFP  np  1  FI  1  T  1  A  |  ANCLF  INPPFMFMTS 

NflM  MPMFF  OF  IN  |  T  |  Al  ~H(lB  \  J  OmT  a  j  Ynpm  FmfKiT  S 

Ain  I  FI  TI1TAI  MOMMPB  IIF  I  M 1  T  I  A I  PAY*; 


Mill  IMF  MIIMHFP  OF  CALCnMP  (.IMPS  NPCFSSAOY  TO  OIITLINF  STBIICTIIBF 
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13  ABSTRACT 


Simulation  of  seismic  rays  for  a  spherical  earth  and  a  flat  earth  has  been  achieved  in  highly  complex 
models.  Travel  times  and  approximate  amplitudes  of  seismic  waves  can  be  found  for  both  two-  and  three- 
dimensional  models  of  portions  of  the  earth.  In  seismology  and  other  disciplines  ray  construction  customarily 
has  been  applied  to  simplified  geometries.  It  has  lx>cn  necessary  to  assume  that  the  seismic  wave  velocity 
distribution  of  the  earth  was  relatively  uniform  and  symmetric. 

Recently,  however,  the  earth  has  been  found  to  be  more  complex  and  non -uniform  than  formerly  assumed. 

A  need  has  thus  arisen  In  seismology  to  test  highly  heterogeneous  models  of  seismic  velocity  distribution. 

At  the  same  time  the  development  of  the  modern  digital  computer  has  provided  a  means  of  performing  the 
necessary  ray  constructions  and  numerical  calculations. 

The  problem  of  complicated  seismic  velocity  distributions  was  therefore  investigated  in  terms  of  the  most 
appropriate  use  of  the  digital  computer.  For  this  investigation  a  velocity  field  was  set  up,  and  the  propagation 
computations  made  for  short  segments  of  rays  within  this  field.  Total  travel  times  arc  found  by  adding  the 
travel  times  oi  connected  ray  segments.  Essentially,  the  nature  of  propagation  was  duplicated  cn  the  com¬ 
puter,  ir.  that,  at  the  location  of  each  segment  along  the  path  of  propagation,  the  initial  condition  and  effect  of 
the  surroundings  determine  the  succeeding  direction  of  the  following  segment. 

Doth  visual  and  numerical  results  have  shown  that  this  simulation  method  can  be  usefully  applied  to  in¬ 
vestigation  of  seismic  velocity  distributions  of  portions  of  the  earth  of  any  size  or  complexity. 
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